In [1]:
using Revise
In [2]:
using PyPlot, LinearAlgebra
pygui(false)
Out[2]:
In [3]:
include("psfeig.jl")
include("utils.jl")
Out[3]:
In [96]:
path = "/home/bruno/Work/Luana_PSF/psfs_su5_size71/all_psfs_su5_size71_com_ph_interp.rsf@"
image = utils.read_bin(path,960,5300);
In [97]:
utils.sisshow(image, vmin=-0.1, vmax=0.1)
In [98]:
b1=70; b2=70; k1=36; k2=36; l1=920; l2=5260;
println("psf grid: ($(length(k1:b1:l1)),$(length(k2:b2:l2)))")
psfs = psfeig.psf_chop(image,b1,b2,k1,k2,l1,l2);
In [99]:
size(psfs.val,3)
Out[99]:
In [100]:
utils.sisshow(psfs.val[:,:,9*13+1], psfs.val[:,:,9*13+5], psfs.val[:,:,9*13+9], psfs.val[:,:,9*13+13], perc=99.9)
In [11]:
o=36; l=12
imshow(psfs.val[(o-l):(o+l),(o-l):(o+l),1], cmap="Greys")
grid()
In [12]:
size(psfs.val[(o-l):(o+l),(o-l):(o+l),1])
Out[12]:
In [13]:
findmax(psfs.val[(o-l):(o+l),(o-l):(o+l),1])
Out[13]:
In [101]:
function plot_psf_with_condition_number(ipsf)
P = psfs.val[(o-l):(o+l),(o-l):(o+l),ipsf];
G = psfeig.psf_to_matrix(n,P);
figure()
imshow(psfs.val[(o-l):(o+l),(o-l):(o+l),ipsf],
interpolation="bicubic", cmap="Greys", vmin=-1, vmax=1)
colorbar()
grid()
title("PSF #$(ipsf): condition number = $(cond(G))")
end
Out[101]:
Distict PSFs
In [15]:
o=36; l=12; n=25;
for ipsf in [9*13+1, 9*13+5, 9*13+9, 9*13+13]
plot_psf_with_condition_number(ipsf)
end
In [17]:
o=36; l=12; n=25;
for ipsf in [9*13+1, 10*13+1, 9*13+2, 10*13+2]
plot_psf_with_condition_number(ipsf)
end
Distinct PSFs
In [102]:
o=36; l=12; n = 100
P1 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+1];
P2 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+5];
P3 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+9];
P4 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+13];
In [103]:
G = psfeig.psf_to_matrix(n,P1,P2,P3,P4); varinfo(r"G")
Out[103]:
In [27]:
cond(G)
Out[27]:
In [28]:
imshow(reshape(G[(75-1)*100+75,:],100,100)); grid()
Neighbour PSFs
In [56]:
o=36; l=12; n = 100
P1 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+1];
P2 = psfs.val[(o-l):(o+l),(o-l):(o+l), 10*13+1];
P3 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+2];
P4 = psfs.val[(o-l):(o+l),(o-l):(o+l), 10*13+2];
In [57]:
G = psfeig.psf_to_matrix(n,P1,P2,P3,P4); varinfo(r"G")
Out[57]:
In [31]:
cond(G)
Out[31]:
In [32]:
imshow(reshape(G[(50-1)*100+50,:],100,100)); grid()
In [39]:
n = 25;
P = psfeig.bluker(12,1.0);
In [107]:
G = psfeig.psf_to_matrix(n,P); varinfo(r"G")
G2 = psfeig.psf_to_matrix(n,P,P,P,P); varinfo(r"G2")
G3 = psfeig.psf_to_matrix(n,P,0.5*P,0.2*P,0.1*P); varinfo(r"G3")
Out[107]:
In [109]:
cond(G)
Out[109]:
In [110]:
cond(G2)
Out[110]:
In [111]:
cond(G3)
Out[111]:
In [105]:
imshow(reshape(G[(11-1)*25+11,:],25,25)); grid()
In [103]:
imshow(reshape(G4[(11-1)*25+11,:],25,25)); grid()
In [40]:
imshow(P)
Out[40]:
In [41]:
findmax(P)
Out[41]:
In [42]:
size(P)
Out[42]:
| $n$ | Time | Size |
|---|---|---|
| 64 | 27s | 128Mb |
| 96 | 162s | 324Mb |
In [36]:
n=100; x = rand(Float32, n*n,n*n); varinfo(r"x")
Out[36]:
In [83]:
@time v = eigvals(x);
In [84]:
x = nothing
In [85]:
GC.gc()
In [35]:
N = 50;
k = 12;
nn = 11;
c = zeros(nn);
sigma = range(0.1, stop=2.0, length=nn);
@time for i = 1:nn
P = psfeig.bluker(k,sigma[i]);
G = psfeig.psf_to_matrix(N,P);
c[i] = cond(G);
end
semilogy(sigma,c);
In [22]:
using FileIO
In [53]:
gulfacks = load(File(format"PNG", "gulfacks.png"));
refl = [ Float32(gulfacks[i,j].r) for i=1:100,j=1:100];
In [43]:
imshow([refl[i,j].r for i=1:100,j=1:100], cmap="Greys")
colorbar()
Out[43]:
In [62]:
function hprod(G,r)
n1,n2 = size(r)
m1,m2 = size(G)
@assert(m1 == n1*n2, "Hessian and reflectivity are not compatible")
return reshape(G*reshape(r,n1*n2),n1,n2)
end
function hinv(G,r)
n1,n2 = size(r)
m1,m2 = size(G)
@assert(m1 == n1*n2, "Hessian and reflectivity are not compatible")
return reshape(G \ reshape(r,n1*n2),n1,n2)
end
Out[62]:
In [69]:
hr = hprod(G,refl);
In [55]:
imshow(hr, cmap="Greys")
Out[55]:
In [105]:
hr = hprod(G,refl);
In [106]:
imshow(hr, cmap="Greys"); title("Com PH")
Out[106]:
In [70]:
refl_inv = hinv(G,hr);
imshow(refl_inv, cmap="Greys"); colorbar()
Out[70]:
In [58]:
hr = hprod(G,refl);
imshow(hr, cmap="Greys")
Out[58]:
In [63]:
refl_inv = hinv(G,hr);
imshow(refl_inv, cmap="Greys"); colorbar();
Out[63]:
In [66]:
Out[66]:
In [72]:
using IterativeSolvers
In [90]:
refl_inv = lsqr(G, hr[:], maxiter=100, verbose=true);
In [107]:
refl_inv = lsqr(G, hr[:], maxiter=100, verbose=true);
In [108]:
imshow(reshape(refl_inv,100,100), cmap="Greys", vmin=0, vmax=1); colorbar(); title("Com PH")
Out[108]:
In [93]:
imshow(reshape(refl_inv,100,100), cmap="Greys", vmin=0, vmax=1); colorbar();
In [ ]: